\documentclass[10pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{url}
\usepackage{framed, color}
\definecolor{shadecolor}{rgb}{1,1,1}
\usepackage{graphicx}
\usepackage{listings}
\usepackage{color}
\usepackage{listings}
\lstset{ %
language=R,                % choose the language of the code
basicstyle=\footnotesize,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\footnotesize,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=7pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{white},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,           % adds a frame around the code
tabsize=2,          % sets default tabsize to 2 spaces
captionpos=b,           % sets the caption-position to bottom
breaklines=true,        % sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={\%*}{*)},          % if you want to add a comment within your code
xleftmargin=.5in,
xrightmargin=.25in
}

\title{Week 2: Single Neuron Models}
\author{Phong Le, Willem Zuidema}

\begin{document}
\lstset{language=R}
\renewcommand{\lstlistingname}{Code}

\maketitle

In today's computer lab we will have a closer look at two differential
equation models of single neuron dynamics: the Fitzhugh-Nagumo model
and the Izhikevich model. We start with a very brief
introduction to dealing with vectors in R and plotting your results,
and then look at how we can use a package for analyzing (linear)
systems of differential equations.

Note: if you are familiar with R you may skip the section~1.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Basic R: vector and graphics}
\label{section basic R}
If you are absolutely new to R, you might want to first go through the
first three pages of this R tutorial:
\begin{itemize}
    \item \url{http://www.illc.uva.nl/LaCo/clas/clc13/assignments/deboer13rtutorial.pdf}
\end{itemize}

Ready? Then, say hello to R by typing \texttt{print("Hello R!")}. 

\subsection{Vector}
There are many ways to construct a vector. The first, also the simplest, way is
to use the \texttt{c} function, e.g., \texttt{vec <- c(1,2,3)} creates a 
\textit{column} vector 
\[
    vec = \begin{pmatrix} 
            1 \\
            2 \\
            3 
        \end{pmatrix}
\]
Another way is to use the \texttt{seq(n,m,by=k)} function, which creates 
$(n,n+k,n+2k,...,n+uk)^T$ vector (where $u=(m-n) \text{ div } k$). For instance, 
the output of \texttt{seq(1,1.5,by=0.1)} is the vector $(1,1.1,1.2,...,1.5)^T$.

\begin{framed}
Exercise \ref{section basic R}.1: Create the following vectors
\[
    a = \begin{pmatrix} 
            1.5 \\
            -1 \\
            3 
        \end{pmatrix} \;\;
    b = \begin{pmatrix} 
            -1 \\
            -1 \\
            -1
        \end{pmatrix} \;\;
    c = \begin{pmatrix} 
            2 \\
            4 \\
            6 \\
            ... \\
            24
        \end{pmatrix} \;\;
\]
\end{framed}

\begin{framed}
Exercise \ref{section basic R}.2: Given vectors in the above exercise, 
\begin{itemize}
    \item What do you get with \texttt{a + b}, \texttt{a - b}, \texttt{a * b},
    \texttt{a / b}?
    
    \item What do you get with \texttt{a > b}, \texttt{a < b}, \texttt{a == b},
    \texttt{a != b}?  
\end{itemize}
\end{framed}

We can get an element by using the operator \texttt{[]}. For instance, 
\texttt{vec[i]} points to the i-th element of a vector \texttt{vec}, 
The operator \texttt{[]} can do more than that: we can get a set of elements.
For instance, \texttt{a[c(1,3)]} and \texttt{a[c(TRUE,FALSE,TRUE)]} 
point to the first and the third elements of vector \texttt{a}.

\begin{framed}
Exercise \ref{section basic R}.3: Find all positive elements in \texttt{a}. 
(Hint: use \texttt{a > 0}.)
\end{framed}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Graphics}

\paragraph{Line plot}
In order to draw a function $y=f(x)$, we can use \texttt{plot(x,y,type='l')} where
\texttt{x} is a vector containing $n$ ordered values of $x$ and \texttt{y} is a
vector such that \texttt{y[i] = f(x[i])}. For instance, the following code will draw
the function $y=sin(x)$
\begin{lstlisting}
x <- seq(-5,5,by=0.1)           # create x = (-5,-4.9,-4.8,...,4.9,5)
y <- sin(x)                     # compute y values
plot(x,y,type='l',col='blue')   # draw the graph, using blue color
savePlot("sin.png",type="png")	# save the current plot to a png file
\end{lstlisting}
If you want to draw another function on the same graph, you need to use 
the \texttt{lines} function. For instance, insert \texttt{lines(x,cos(x))}
before the last line. Now, you should get a graph similar to 
Figure~\ref{fig sin cos}.
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.5\textwidth]{sin.png}
    \caption{Drawing two functions, $y=sin(x)$ and $y=cos(x)$ in one graph.}
    \label{fig sin cos}
\end{figure}


\paragraph{Scatter plot}
Remove \texttt{type='l'} and replace \texttt{lines(x,cos(x))} by
\texttt{points(x,cos(x))}, what do you get?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ODE and Phase Plane Analysis}
\label{section ode}

\paragraph{Required Libraries}
There are two libraries we need for the task, \texttt{deSolve} for
solving ODEs, and \texttt{fields} for drawing vector fields. Install them by 
typing
\begin{lstlisting}
install.packages("deSolve")
install.packages("fields")
\end{lstlisting}
If you install a library on a computer where you don't have administrator 
rights, e.g., a university computer, you need to use
\texttt{install.packages("library\_name", lib="your\_lib\_path")}
where `your\_lib\_path' is any directory where you have rights to write files.
In order to use a library, execute \texttt{library(lib\_name)}, for instance 
execute \texttt{library(deSolve)}.

\paragraph{Required R Code} At \url{http://www.illc.uva.nl/LaCo/clas/fncm13/assignments/computerlab-week2/} you can find the four R-files you need for this exercise: 

\begin{itemize}
	\item `draw.R' contains functions for drawing phase planes and
          trajectories,
	
	\item `linear\_ode.R', `FN\_ode.R' and `I\_ode.R' contain R code 
	for linear models, the Fitzhugh-Nagumo model, and the Izhikevich model.
\end{itemize}
You will be asked to add your own R code to those files.

\subsection{ODE}
A system of first-order ODEs we are interested in is defined as follows
\begin{align*}
    \frac{du_1}{dt} &= f_1(t,u_1,...,u_n,p_1,...,p_m) \\
    \frac{du_2}{dt} &= f_2(t,u_1,...,u_n,p_1,...,p_m) \\
    ...&& \\
    \frac{du_n}{dt} &= f_n(t,u_1,...,u_n,p_1,...,p_m)
\end{align*}
where $t$ is an independent variable (e.g., time), $u_1,...,u_n$ are
dependent variables (e.g., variables describing the state of a dynamic
system), and $p_1,...,p_m$ are parameters. 

%Now, before turning to the next part, we strongly recommend you to read
%through Chapter 11, ``Introduction to Phase Plane Analysis'', in the book
%Wallisch et al. 
%``MATLAB for Neuroscientists: An Introduction to Scientific Computing in\
%MATLAB'' in order to get many important 
%concepts such as \textit{phase plan}, \textit{equilibrium}, \textit{fix point}, 
%\textit{nullcline}, \textit{vector field}, \textit{trajectory}, etc.

\subsection{Phase Plane Anaylysis}
In this section, we will consider a system of two linear first-order ODEs
\begin{align*}
    \frac{dx}{dt} &= ax + by \\
    \frac{dy}{dt} &= cx + dy
\end{align*}
where $x,y$ are dependent variables, and $a, b, c, d$ are parameters. 
This system has two linear nullclines:
\begin{align*}
    \frac{dx}{dt} &= ax + by = 0 \\
    \frac{dy}{dt} &= cx + dy = 0
\end{align*}
The intersection of the two nullclines, if it exists, is called \textit{equilibrium}.

Now, open the file `linear\_ode.R', you should find the following function
which declares the system of ODEs above:
%% # a system of two linear first-order ODEs
%% # dx/dt = ax + by
%% # dy/dt = cx + dy
\begin{lstlisting}
# input:
#   t: current time point
#   var: c(x,y) current state
#   p: params, c(a,b,c,d)
# output: list(c(dx/dt,dy/dt))
linear.ode.system <- function(t, var, p) {
    x <- var[1]; y <- var[2]
    a <- p[1]; b <- p[2]; c <- p[3]; d <- p[4]

    return(list(c(a*x + b*y, c*x + d*y)))
}
\end{lstlisting}
and the following function for nullclines
\begin{lstlisting}
# give x, find y such that
# dx/dt = dy/dt = 0
linear.nullcline <- function(x, p) {
    a <- p[1]; b <- p[2]; c <- p[3]; d <- p[4]

    return(c(-a/b*x, -c/d*x))
}
\end{lstlisting}

\begin{framed}
Exercise \ref{section ode}.1: Check for yourself how these functions work.
\end{framed}


%% \begin{lstlisting}
%% # draw phase plane
%% # input:
%% #   ode.sys:  ode system
%% #   ode.event, ode.root: event and root functions for deSolve 
%% #                       in the case that ode.sys is discontinuous
%% #   nullcline: function to compute nullclines
%% #   params: params for ode system
%% #   xrange: c(xmin,xmax) of the graph
%% #   yrange: c(ymin,ymax) of the graph
%% #   xstep, ystep): the x, y distance between two vector tails
%% #   x.ini, y.ini: initial conditions of the system
%% #   times: time sequence of trajectories
%% draw.phase.plane <- function(ode.sys, ode.event, ode.root,
%%                     nullcline, params,
%%                     xrange, yrange, xstep, ystep,
%%                     xlabel, ylabel,
%%                     x.ini, y.ini, times)

%% # draw trajectory over time
%% # input:
%% #   ode.sys: system of ODEs
%% #   params: parameters for ode.sys
%% #   x.ini, y.ini: initial condition
%% #   times: time steps
%% #   xlabel, ylabel: names of the first and the second dependent variables
%% draw.trajectory.over.time <- function(ode.sys, ode.root, ode.event, params,
%%                                         x.ini, y.ini, times,
%%                                         xlabel, ylabel)            
%% \end{lstlisting}

Now also look at the functions draw.phase.plane and
draw.trajectory.over.time in the file `draw.R'. Note that you don't
need to fully understand the bodies of the functions, just skim them
to get an idea how they work. Execute the code in draw.R (by using the
menu or by typing \verb+source(file.choose())+ in the console and
selecting draw.R).

To get an idea what the arguments are, let's execute
`linear\_ode.R'. To do that, change your working directory in the menu
to the folder where linear\_ode.R is and then execute it with
\verb+source("linear_ode.R")+.

\begin{lstlisting}
######### linear ode ###########
test.ode.system = linear.ode.system
test.nullcline  = linear.nullcline
test.params     = c(4,-3,15,-8)
test.xrange     = c(-10,10)
test.yrange     = c(-10,10)
test.xstep      = 1
test.ystep      = 1
test.xlabel     = "x"
test.ylabel     = "y"
test.x.ini      = c(2)
test.y.ini      = c(-10)
test.times      = seq(1,5,by=0.01)
test.ode.root   = NULL
test.ode.event  = NULL

########### run experiment ##########
draw.phase.plane(
            ode.sys = test.ode.system, 
            ode.event = test.ode.event, ode.root = test.ode.root,
            nullcline = test.nullcline,
            params = test.params,
            xrange = test.xrange, yrange = test.yrange,
            xstep = test.xstep, ystep = test.ystep,
            xlabel = test.xlabel, ylabel = test.ylabel,
            x.ini = test.x.ini, y.ini = test.y.ini, times = test.times)

draw.trajectory.over.time(
                        ode.sys = test.ode.system, 
                        ode.event = test.ode.event, ode.root = test.ode.root,
                        params = test.params,
                        x.ini = test.x.ini, y.ini = test.y.ini,
                        times = test.times,
                        xlabel = test.xlabel, ylabel = test.ylabel)
\end{lstlisting}
The first four arguments (lines 19-21) are to tell which system of ODEs we are interested in. 
The fifth argument (line 22), \texttt{params = c(4,-3,15,-8)}, 
says that the parameters for the system are
\texttt{a = 4, b = -3, c = 15, d = -8}. Some of the other arguments are described in
Figure~\ref{fig linear odes}.
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.7\textwidth]{linear.png}
    \caption{Phase plane of the system of two linear first-order ODEs with 
    \texttt{a = 4, b = -3, c = 15, d = -8}.}
    
    \label{fig linear odes}
\end{figure}

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.5\textwidth]{trajectory.png}
    \caption{$x$ and $y$ versus $t$.}
    
    \label{fig trajectory}
\end{figure}

\begin{framed}
Exercise \ref{section ode}.2: With the system given in Figure~\ref{fig linear odes}
and \ref{fig trajectory},
\begin{itemize}
    \item Show that (0,0) is an equilibrium.
    \item Explain why, if the initial condition is a point very close to (0,0),
    the system will definitely end up at that point.
\end{itemize}
\end{framed}

\begin{framed}
Exercise \ref{section ode}.3: Modify the file `linear\_ode.R' to 
draw the phase plane of the system
\begin{align*}
    \frac{dx}{dt} &= x + y \\
    \frac{dy}{dt} &= 4x + y
\end{align*}
with the two initial conditions (4.8,-10) and (5.1, -10). 
(Tip: you can plot multiple trajectories in the same graph by letting text.x.ini and text.y.ini contain a list of x- and y-coordinates. E.g., 
\begin{verbatim}
	test.x.ini		= c(-3,3,1.2)
	test.y.ini		= c(3,-3,-0.9)
\end{verbatim})
\begin{itemize}
    \item Explain why the two trajectories have those shapes.
    \item Find an initial condition such that the trajectory ends up at the equilibirum. 
\end{itemize}
\end{framed}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Fitzhugh-Nagumo Model}
\label{section FN}
In this section, we will explore the Fitzhugh-Nagumo model, which is described 
by the following system of two ODEs
\begin{align*}
    \frac{dv}{dt} &= c(v-\frac{1}{3}v^3 + r + I) \\
    \frac{dr}{dt} &= -\frac{1}{c}(v - a + br)
\end{align*}

\begin{framed}
  Exercise \ref{section FN}.1\footnotemark: (Use the file `FN\_ode.R';
  to change the parameters, you will need to change the line
  \verb+test.params = c(0.7,0.8,3,0)+. Tip: put a \verb+#+ in front of
  test.params in FN\_ode.R. You can then type
  \verb+test.params = c(0.7,0.8,3,-0.2)+ etc. in the console, and run
  the script with \verb+source("FN_ode.R")+.)

\begin{itemize}
    \item Draw the phase plane of the Fitzhugh-Nagumo model with 
    $a=0.7,b=0.8,c=3,I=0$. Is the equilibrium point stable (i.e., 
    are trajectories attracted to this point or repelled from it)?
  \item Try out different initial values for the trajectories (by changing text.x.ini and text.y.ini). 
  \item Try 'depolarizing' the neuron from its equilibrium value at $(v, r) = (1.1994, -0.62426)$. What happens to the trajectories starting from $v=1, 0.5, 0.0, -0.5$? I.e., try:
\begin{verbatim}
	test.x.ini		= c(1, 0.5, 0.0, -0.5)
	test.y.ini		= c(-0.62,-0.62,-0.62,-0.62)
\end{verbatim}

    \item Change the injected current value to $I=-0.2$, the initial condition 
    is $(v, r) = (1.1994, -0.62426)$. Is this point still stable?
    
    \item Determine what $v$ versus $t$ looks like for a trajectory on this 
    phase plane. Would you classify the injected input of -0.2 as a 
    superthreshold or subthreshold stimulus? Does this neuron exhibit 
    subthreshold oscillations for this value of injected current?
    
    \item Change the injected current value to $I=-0.4$ with the initial 
    condition $(v, r) = (1.1994, -0.62426)$. Is this point still stable? 
    Plot several trajectories on this phase plane. Since the nullclines 
    intersect at only a single point, there are no other equilibrium points 
    for this system, but trajectories may be attracted to some other closed 
    orbit: a limit cycle. Is there such a limit cycle in this system? 
    
    \item Finally, repeat the analysis for $I=-1.6$ and examine $v$ versus $t$.
    Does this neuron spike continuously as it did before? Neurons are known to 
    exhibit a phenomenon called excitation block, whereby increasing the current 
    injection can often repress repetitive firing behavior.
\end{itemize}
\end{framed}
\footnotetext{This exercise is the project in 
Chapter 12, ``Exploring the Fitzhugh-Nagumo model'', in the book
Wallisch et al.
``MATLAB for Neuroscientists: An Introduction to Scientific Computing in\
MATLAB''.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Izhikevich Model}
\label{section I}

In this section, we will explore the Izhikevich model, which is described 
by the following system of two ODEs
\begin{align*}
    \frac{dv}{dt} &= 0.04v^2 + 5v + 140 - u + I \\
    \frac{du}{dt} &= a(bv - u)
\end{align*}
with the reset condition: $\text{if } v \ge 30 \text{ then } 
v \leftarrow c;\;\; u \leftarrow u+d$. 
Because of the reset condition, the model turns out discontinuous at the 
reset point, and modelling it as an input to \texttt{deSolve} is a bit tricky.
Therefore, in this section, we give you Izhikevich model code, and your
task is simply try out some parameter sets.

\begin{framed}
  Exercise \ref{section I}.1: (Use the file `I\_ode.R')
  Plot the phase space of the model with parameters $(a, b, c, d, I) =
  (0.02, 0.2, -65, 8, 0)$, with null clines and a couple of
  trajectories that start from various points around the
  equilibrium. Do you observe qualitatively similar behavior to the
  Fitzhugh-Nagumo model with $I=0$?
\end{framed}

\begin{framed}
Exercise \ref{section I}.2\footnotemark: (Use the file `I\_ode.R')
The purpose of this exercise is to see if you can discover what parameter 
sets lead to regular spiking, fast spiking, or intrinsically bursting behavior. 
What kinds of behaviors do the following parameter sets produce? (with $I=10$)
\begin{itemize}
    \item $(a, b, c, d, I) = (0.02, 0.2, -65, 8, 10)$
    \item $(a, b, c, d, I) = (0.02, 0.2, -55, 4, 10)$
    \item $(a, b, c, d, I) = (0.1, 0.2, -65, 2, 10)$
    \item $(a, b, c, d, I) = (0.1, 0.25, -65, 2, 10)$
\end{itemize}
(You may need to check this 
\url{http://www.izhikevich.org/publications/spikes.htm}.)
\end{framed}
\footnotetext{This exercise is Exercise 22.1 in 
Chapter 22, ``Simplified Model of Spiking Neurons'', in the book
Wallisch et al.
``MATLAB for Neuroscientists: An Introduction to Scientific Computing in\
MATLAB''.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Submission}

You have to submit a file named `your\_name.pdf' for those exercises
requiring explanations and math solutions through Blackboard before
15:00 Monday 11 Nov.

\end{document}
